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Abstract. Condensation transition in a non-Markovian zero-range process is stud- 
ied in one and higher dimensions. In the mean-field approximation, corresponding to 
infinite range hopping, the model exhibits condensation with a stationary condensate, 
as in the Markovian case, but with a modified phase diagram. In the case of nearest- 
neighbor hopping, the condensate is found to drift by a "slinky" motion from one site 
£H ■ to the next. The mechanism of the drift is explored numerically in detail. A modified 

model with nearest-neighbor hopping which allows exact calculation of the steady state 
is introduced. The steady state of this model is found to be a product measure, and 
the condensate is stationary. 
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1. Introduction 

\ In recent years, much progress has been made in the theoretical understanding of 



nonequilibrium condensation |lH3[. Condensation phenomena of this type are known to 
occur in a large variety of systems, where a macroscopic fraction of a conserved "mass" 
is accumulated in a microscopic portion of an extended system. Some examples of such 
systems are compartmentalized shaken granular gases |4], gelation, i.e., the formation 



of a macroscopic hub in complex network |5j, condensation of wealth in economics 
and the formation of traffic jams on highways [3]. 

An elucidation of the mechanism which lies behind the condensation transition in 
such systems was achieved by studying simplified but prototypical toy models, most 
notably the zero-range process (ZRP). In this process, particles hop stochastically 
between boxes with hopping rates which depend only on the occupation of the box 
from which each particle departs. This models diffusing particles which "interact" only 
with particles at the same location (i.e., the same box), and hence the name of the 
process. For any choice of hopping rates, the steady-state distribution of particles is 
known to factorize into single box terms and thus it may be computed exactly. Using 



Motion of condensates in non-Markovian zero-range dynamics 



2 



the exact solution, it was shown in [8] that a condensation transition may occur when 
the rate of hopping out of a box decreases with its occupation (modeling an attractive 
"interaction"). In this condensation transition, the system is homogeneous at small 
particle densities, while when the density exceeds a critical value, a single (randomly 
chosen) site is occupied by a macroscopic fraction of all particles. 

When studying systems which exhibit a condensation transition, such as those 
listed above, the ZRP may be used to gain qualitative, and sometimes also quantitative, 
insight. This is done by mapping the dynamics of the system under consideration, 
usually in a n ap proximate way, to that of the ZRP, and then utilizing known results for 
the ZRP 0-E3- Such a mapping is achieved by disregarding some of the structure of 
the original system, so that it can be reduces to the simplified balls and boxes picture of 
the ZRP. Underlying this procedure is an assumption that the condensation transition 
in the ZRP is universal in some sense, i.e., that the details which are lost when mapping 
a system to the ZRP are irrelevant to condensation. 

This universality assumption was recently examined by studying how disregarding 



dynamical degrees of freedom may affect the condensation transition [14J. There, a small 
variation of the ZRP, which introduces temporal correlations in its dynamics, was found 
to have two main effects on the condensation transition: (i) a calculation in a mean- 
field setting showed that the non-Markovian nature of the dynamics "renormalizes" the 
parameters and critical exponents of the ZRP, and (ii) a numerical study of the model 
on a one-dimensional ring revealed that the temporally-correlated dynamics induces a 
"slinky" motion of the condensate throughout the system, whereby the condensate spills 
over from one site to the next. Thus, the nature of the condensed phase is modified in 
a qualitative manner. 

In this paper we present a detailed analysis of the model introduced in [14|. First, 
we elaborate on the mean-field solution which was presented in 14] and generalize the 
results to a much broader class of hopping rates. We then turn to dynamics on a ring 
with hopping to the nearest-neighbor site and present a detailed study of the mechanism 
for the condensate drift. In studying the behavior of finite rings (of size L < 2000 sites), 
we identify two different modes of condensate motion: motion through a barrier, in 
which the condensate is carried by a single site and spilling to the next site is initiated 
only once it overcomes a barrier, and motion with no barrier, in which the condensate 
is in a continual motion. 

In addition, we consider a somewhat modified non-Markovian dynamics which 
allows an exact computation of the steady state on lattices of any dimension, even 
in the case of nearest-neighbor hopping. For this variant of the model we show that 
the steady-state distribution is given by a product measure, like in the Markovian ZRP, 
even though the currents of particles are temporally correlated. The exact solution 
of the model is shown to be qualitatively similar to the mean-field calculation, and in 
particular there is no condensate drift. 

The paper is organized as follows. After describing the non-Markovian ZRP in 
Sec. [2J we present in Sec. |3] the full mean-field solution of the model and study how 
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condensation is affected by the non-Markovian dynamics. The numerical study of the 
condensation transition on a lattice with nearest-neighbor hopping is presented in Sec. 
HI There we examine a particular choice of rates, which we term the on-off model, on a 
symmetric, totally asymmetric and partially asymmetric dynamics, and we analyze the 
mechanism for the motion of the condensate. In Sec. El we present the exactly solvable 
variant of the non-Markovian ZRP, calculate its steady state distribution and discuss 
condensation in this model. 

2. Description of the model 

We consider a system of N particles hopping between L boxes (labeled by i — 1, . . . , L) 
with a mean density p = N/L. We are mainly interested in the thermodynamic limit 
in which L, N — > oo with the density p kept constant. The state of each box is given 
by two variables: the number of particles in the box n«, and a "clock" variable 7$. Both 
variable take non-negative integer values: n^, = 0, 1, 2, . . .. A configuration of the 
system is thus given by the set of pairs (n, r) = Ti)}f =1 . 

The dynamics proceeds by particles jumping one at a time between the boxes, and, 
in parallel, by advances of the clocks. Conforming with the "zero-range" character of 
the ZRP, the rate with which a particle leaves a box is taken to depend only on the 
state of the box, i.e., on its occupation and clock state. We denote these hopping rates 
by u(n, r). The clock dynamics is correlated with particle jumps: every time a particle 
jumps into a box, its clock is reset to zero. Independently, the clocks are advanced with 
a constant probability per unit time c. The two types of dynamical moves (a jump of a 
particle between two boxes % and j, and and an advance of the clock at site i) can be 
written as 

(n h n), (rij, Tj) PlJ !^' Tl) (nj-1, n), (nj + 1, tj = 0) 

(rii,Ti) (ni,Ti+i). (1) 

Here, Pij IS cL connectivity matrix which states the probability that a particle departing 
from site i will choose site j as its target. Particular choices of Py are further discussed 
below. The jump rates must satisfy u(0, r) = for all r, as a particle cannot jump out 
of an empty box. 

To be precise about the meaning of "rate" this dynamics can be rephrased as follows: 
each box % of the lattice carries two alarm-clocks which ring after some random time. 
All 2L alarm-clocks ring independently. Given that the state at box i is (nj, Tj), alarm- 
clock number 1 of this site rings after an exponentially distributed random time with 
parameter u(n,r), while alarm-clock number 2 rings after an exponentially distributed 
random time with parameter c. If clock 1 rings first, a particle selects a target box j 
with probability and jumps to it. If clock 2 rings first, the internal clock at box % is 
incremented by one unit. After any change of the state at box i, the clocks ring again 
after an exponentially distributed random time determined by the updated state at box 
i. 



Motion of condensates in non-Markovian zero-range dynamics 



4 



The dynamical rules ([I]) are written for a general connectivity matrix p^. In 
this paper we concentrate mainly on two schemes for the choice of target box: mean- 
field (MF) dynamics, and a one- dimensional homogeneous ring geometry with nearest- 
neighbor hopping. In mean-field (MF) dynamics, the target box j ' ^ % is chosen randomly 
and uniformly from all boxes, i.e., 

Pij = J~[ for a11 ^ *■ ( 2 ) 

In the case of ring dynamics, on the other hand, particles hop only between nearest- 
neighbor sites, possibly in an asymmetric fashion. Accordingly, the target box is chosen 
to be % + 1 with probability 1 — p and % — 1 with probability p (where box L + 1 is 
identified with box 1), i.e., 

r (i-p) ^ j =2+i) 

Pij = I P (if j = i ~ 1) ■ (3) 

I otherwise 

Here, < p < 1 is the asymmetry parameter: when it is equal to the hopping is totally 
asymmetric and no hopping backwards can occur, while when p = 1/2 the dynamics is 
completely symmetric. Below, MF dynamics is studied in Sec. EJ while ring dynamics is 
studied in Sections EH and [5j We also briefly examine (in Sections l4.2.3l and 15. 2p dynamics 
on higher- dimensional lattices. The generalization of fl3]) to the higher-dimensional case 
is rather straightforward and will be presented below when it is discussed. 

The dependence of the jump rates on r renders the particle jump process, when 
it is taken by itself, non-Markovian, as the rate of a jump depends on how much time 
has passed since a particle hopped into the jump site. When considered in the higher 
dimensional space of occupations together with clocks, the full jump /increment process 
defined above is Markovian. Nevertheless, we will refer to this process as the non- 
Markovian ZRP, to stress the history dependence of the jump process. 

The process may be implemented by a discrete-time Monte-Carlo version of 
this dynamics with random sequential update which is defined as follows: define 
Pmax = max„ iT [M(n, r) + c]. For the Monte-Carlo update pick a random box uniformly 
and attempt to make one of the following changes: (i) move a particle to a target 
box (selected according to the appropriate scheme) with probability u(n,r)/p max , (ii) 
increment the internal clock with probability c/p max . A total of Lp m3iX consecutive 
update attempts constitute one Monte-Carlo time unit. 

A final remark on the nature of the clock variables. As is clear from the dynamical 
rules, the clock variables proceed in an irregular, stochastic fashion. Therefore, they 
do not measure the exact time that has passed since a particle last entered each site. 
Choosing clock variables which really measure time, i.e., which are continuous and 
proceed regularly, might seem more natural for some physical applications. Such regular 
clocks are not considered below, but we remark that they may be achieved starting 
from the dynamics (pQ) by taking an appropriate limit. This procedure is described in 
Appendix A 



Motion of condensates in non-Markovian zero-range dynamics 



5 



3. The non-Markovian ZRP with mean-field dynamics 

3.1. General observations 

In this section we analyze the non-Markovian ZRP with mean-field dynamics. In 
particular, we investigate how condensation is affected by the non-Markovian nature 
of the jump process. 

Since each jump of a particle is to an arbitrarily chosen box, the MF dynamics does 
not generate correlations between different boxes beyond the correlations which arise 
from conservation of particles. Therefore, in the thermodynamic limit, the stationary 
distribution is expected to factorize into a product of single-box terms 

L L 

V(n, t) = Z L ] N J] P(rk, n) 6 m - N) , (4) 

i=l i=l 

where P(ni, 7$) are the single-box occupation and clock probabilities, and the 6 function 
is a consequence of the conservation of particles. The normalization is given by 

L L 
n,T i=l i=l 

As with the Marokovian ZRP, the factorized stationary distribution provides the 
means for an analytic treatment of the model. In the thermodynamic limit, the single- 
box probabilities in the steady state P(n, r) are equal to those of a single box with 
a "mean-field" incoming current J which is generated by all other sites. The master- 
equation for the single site-box probability is 

dP(n t) 

^-1 = JP(n-l)5 Ti0 + cP(n, T-l)+u(n+l, r)P(n+l, r)-P(n, r)[J + c + u(n, r)}, (6) 

where the marginal occupation distribution is defined by P(n) = P(n, r). The first 
term on the RHS corresponds to the box reaching the state (n, r = 0) by a particle 
entering a box with n — 1 particles, the second to an advance of the clock into state r, 
the third to a particle leaving a box with n + 1 particles, and the last to these three 
processes occurring when the box is in state (n, r). This equation is also valid for n = 
or t = if one defines P(— l,r) = P(n, — 1) = (and, as stated above, u(0, r) = 
must also hold). Once Eq. (|6]) is solved for a given MF current J, the current and the 
probability distribution are obtained by the self-consistency requirement 

J = J2^r)P(n,r). (7) 

In the steady state, dP(n, r)/dt = and the master equation Q yields 

P(n,r)[J+c+u(n,r)} = JP(n-l)S rfi +cP(n,T-l)+u(n+l,T)P(n+l,T).(8) 

Summing over all values of r, the terms containing c drop out telescopically, and one is 
left with the recursion relation 

u{n)P(n) = JP(n-l) (9) 
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Here u(n) is the mean hopping rate out of a site with n particles 

Equation (jHJ) expresses the balance between the probabilities to jump into and out of 
a box with n particles. Iterating relations fl9]) yields, as in the Markovian ZRP, the 
steady-state occupation probability 

P(n) = P(0)r/», (11) 
where the single-site weights are given by 

n 

/»=n^)^ i ) (12) 

fc=i 

and P(O)^ 1 = 1 + ^^Lx J n f(n) ensures the proper normalization of P(n). The marginal 
distribution fTTTj) and (fT2]) is the same distribution one obtains for a Markovian ZRP 
but with the jump rates u(n) replaced by the effective rate u(n) [ll. 

Since the stationary distribution of our model has a similar form to that of a 
Markovian ZRP, the analysis of condensation in the model may also proceeds in a similar 
fashion. We therefore briefly review how condensation takes place in the Markovian 
ZRP [1]. The occurrence of condensation in the Markovian ZRP is determined by the 
asymptotic behavior of the jump rates u{n) for large n. Two types of condensation may 
be distinguished: strong condensation, which occurs when the hopping rates tend to zero 
for large n, and weak condensation, which may occur when the hopping rates decrease 
to a constant value. In condensation of the strong type all particles accumulate in one 
box and the current vanishes in the thermodynamic limit. This condensation occurs at 
all densities (i.e., the critical density for condensation is p c = 0). Weak condensation 
takes place only when the rates decrease to a constant more slowly than 1 + 2/n. In 
particular, when the rates have the form 

u ( n ) = 7 (l + b/n a + o(l/n CT )) (13) 

for large n, condensation occurs above some critical density provided that a < 1 or 
a = 1 and b > 2. The critical density p c is non-universal, i.e., it depends on the exact 
form of the rates u(n). Importantly, the parameter 7 only sets the time scale for the 
process and has no affect on the stationary distribution and the condensation transition. 

In the weak condensation scenario, a single site (the condensate), chosen 
spontaneously at random, accommodates L(p — p c ) particles, while the density at all 
other sites remains p c . The condensation transition is thus manifest in the occupation 
probability of a single site, P(n). For the marginal case of a = 1, the probability to find 
n particles in a given site decays exponentially as P{n) ~ n~ b e~ n ^ for densities below 
the critical density, where £(p) diverges as p c is approached. At the critical density, the 
occupation probability has a power law tail P{n) ~ n~ b . Above the critical density, the 
occupation of all background sites remains power-law distributed, while the occupation 
of the condensate is narrowly distributed around L(p — p c ) 16 . 
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3.2. Condensation in an "on-off" model 

As the effective jump rates u(n) play the role of u(n) in the non-Markovian ZRP, it is 
their asymptotic behavior for n — > oo which determines condensation in the model. The 
remainder of this section concentrates on the determination of this asymptotic behavior. 
We begin by discussing a simple choice of jump rates — an "on-off" model which will 
now be introduced — before turning to an analysis of more general jump rates. 
We start the discussion by considering jump rates of the form 

, , f r = ("off" state) 
uM = \u{n) r>l ("on" state). (14) 

In this case, every time a particle hops into a box, that box is turned "off". When 
the box is in this off state no particle can leave it. After an exponentially distributed 
random time (with parameter c) the box is turned back "on", and particles can once 
more jump out of it with a rate u(n). The model with these special rates will be called 
the on-off model. 

In the on-off model, the dynamics depends only on whether r = or r > 1, and 
thus the clock has effectively only two states. Correspondingly, the state of a box can 
be characterized by P s(n) = P(n,r = 0) and P on (n) = J2 T>1 P(n,r). The stationary 
master equation (jSJ) is then given by 

P oS (n)[J + c] =JP(n-l) (15) 
P on {n)[J + u{n)} = cP oS {n) + P on (n+l)u(n+l). (16) 

The solution of these equations is made simple, compared with a general non-Markovian 
ZRP, because the term p(n + l)u(n + 1, 0) which should appear in the RHS of ( TT51) (see 
Eq. (jS])) vanishes. 

To solve these equations we first note that by summing Eq. ( TT5l) over all values of 
n we find that the probability to find a site in the off state is 

P off = J]p off ( n ) = -L (17) 

n 

Next, an expression for P g(n) is found from Eq. ([15]) together with (JSJ) and (fT7|) . 

, , u(n]P(n] P n fr , . . , , . 

Pos(n) = V, = -fu(n)P(n). 18 
J + c J 

A similar expression for P on (n) is found by substituting the rates u(n, r) (which are of 
the form (fl4l) ) into Eq. (TTO]) . yielding 

„ . . u(n) „. . . 
u[n) 

The effective jump rates u(n) can now be obtained from (TIBl) and (Tl9l) using 
Pos( n ) + Pon{n) — P( n ), an d are given by 

1 Pnff 1 

u[n) J u(n) 
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This equation states that the mean time between hops from a site with n particles is 
equal to the mean time this site is in an "off" state plus the time it takes a particle to 
hop out once the system is already "on" . 

Using Eqs. ( TTTj) . ( !T2|) . ( TT7|) and (|20|) . it is now possible to obtain P(n) for any J, 
and subsequently the entire probability distribution is found via ( 1T81) and ( 1T91) . Note 
that P{n) which is found this way depends on J both directly, as seen in Eq. (fTTl) . and 
indirectly through the effective rates (1201) . To finish the calculation, one must find the 
dependence of the current J on the density p. This can in principle be achieved by 
inverting the relation p(J) = ^2 n nP(n). The effective hopping rates ( 120]) are thus a 
function of the density. 

To determine whether or not condensation may occur in the model, only the 
asymptotic form of u(n) is needed. By examining Eq. ( |20l) it is seen that u(n) decreases 
to zero when n — > oo if and only if u(n) decreases to zero, and similarly u{n) decreases 
to a constant if and only if u(n) decreases to a constant. Therefore, strong condensation 
is not affected by the clock-dependent dynamics. To study weak condensation, assume 
jump rates of the asymptotic form (|T3|) with 7 = 1 (as explained above, 7 sets the time 
scale of the process, and can be set to 1 without loss of generality). From Eqs. ( IT71) and 
( 120]) we find, to leading order in 1/n 

^H-^rfl + ^ + of^l (2D 



c + J + 1 V n a \n 
which is again of the form (fT3j) but with an effective hopping parameter 

*- = ttjtt" < "■ < 22) 

If a < 1 condensation occurs in the on-off model as it does in the Markovian ZRP. 
In the commonly encountered case of a = 1, however, condensation only occurs when 
b c s > 2. The critical current at the condensation transition is given in this case by 
J c = lim^oo u(n) [H, which yields, according to (l2~Tj) . J c = (c + J c )/(c + J c + 1), or 



Jc = ~(\l + --1). (23) 



2 VV c 
This allows us to write 



'Jeff 



J c b. (24) 



Therefore, condensation takes place when the hopping parameter satisfies b > 
I + 4/c — . The critical value of b is larger than 2, in contrast with the 
Markovian case for which the critical value for condensation is b = 2. 

3.3. Condensation in MF models with more general rates 

In the previous section we have seen that in the case of jump rates with an asymptotic 
form 

u (n) = 1 + b/n + ..., (25) 
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the on-off dynamics leads to an effective value of b, and thus it may affect the occurrence 
of the condensation transition. We now demonstrate that this holds also when the clock 
dependence is more general than the on-off case, and we show how b e s may be calculated. 
To this end we consider rates of the form 



u(n, t) = u{n)v{r) 



1 + - + ...)v(t) 

n 



(26) 



where u{n) has been taken to be of the form (j25|) . 

As mentioned above, the stationary master equation ([8]) is harder to analyze when 
the rates are not of the on-off type, because P(n, r) depends in ([8]) on P(n,r + 1). 
However, since only the large n asymptotics of u(n) and P(n,r) at criticality affect 
condensation, it is possible to make progress by restricting the discussion to these 
quantities. We therefore assume that J = J c and make the following ansatz: 



u(n) = J c 
P(n,r) = 



hs ( 1 
1 + — + o - 

n \nJ J 



An 



air) 



d(r) (\ 
n \n 



(27) 



P{n) = An 



, d f l 
1 + - + - 

n \n 



This ansatz is motivated by the solution of the on-off model (compare with Eqs. 
( I18p and fl2T]) ). The constant A is a normalization constant, and from the definition 

P{n) = J2 T P( n > T ) it i s seen t na t ^2 T a ( T ) — 1 an d J2 T a ( T )d( T ) = d must hold. 

Substituting the ansatz (1271) in the stationary master equation ([8]) and equating 
terms order by order in 1/n yields to order 0(1) 

J r / c 



a r 



J c + c\J c + c 



and to order 0{l/n) 

d{r) = d- 



J c + c 



i 



5cff- 



(2f 



(29) 



The current and 6 e g can now be found by substituting fl27j) in the definition of u{n) (Eq. 
( TTUl) ) and equating once again order by order in 1/n. To order 0(1), an equation for 
the critical current is obtained 

J r 



E 



T = 



J c + c ^\J C + c 



v(t), 



where fl28|) was used in the last equality. To order 0(l/n), using fl28l and (T59 
found to satisfy 

OO T _ ^ 



(30) 
b e s is 



r=0 

oo 



T=0 



Jc + C 



v{t)v{t') 



(31) 
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Figure 1: The mean- field value of b e s/b as a function of Vq for the rates ( [32]) . The 
line corresponds to the prediction (I3TI) . while the symbols were obtained by numerically 
integrating the mean- field master equation (jfjj) with J = J c of ( J30l) . Note that b e s may 
be either smaller or larger than b, depending on Vq. 

The calculation outlined above is valid as long as the series in ( )30l) and ( 13TT) converge. 
The exponential form of ol(t) in Eq. fl28|) implies that convergence is guaranteed if f (r) 
decays or grows slower than exponentially. In particular, this implies that the results 
are correct if v(t) tends to a finite (non-zero) constant for large r. Note that if v(t) 
decays to zero fast enough, although the series converge any system of a finite size 
will eventually be frozen in an absorbing state in which all r's tend to infinity and no 
particles jump. 

Equation ( 13T1) implies that, as found in the particular case of the on-off model, 
the condensation behavior depends on the memory effects induced by the clocks. Note, 
however, that unlike the on-off case, b e g is not necessarily smaller than b. For instance, 
consider rates of the form (1261) with 



For Vq — these rates reduce to the on-off model, while v o = 1 is the Markovian ZRP. 
For arbitrary t> , Eq. fl30l) yields J c = (v — c + a/ (v — c) 2 + 4c)/2, and Eq. fl3Tj) yields 
b c s = b (c+voJ c )/(c+voJ c + J c — v ). This result, which is plotted in Fig.CEJ demonstrates 
that for different values of v o and c, the effective hopping parameter b e s might be larger 
or smaller than the "bare" value b. 

4. The non-Markovian ZRP with nearest-neighbor dynamics 

The results of the previous section demonstrate that temporal correlations in the 
dynamics of a mean-field ZRP affect the condensation transition. In this section we 
examine whether this mean-field picture persists also when the dynamics allows only 
nearest-neighbor hopping, and whether new effects appear in the latter case. 




(32) 
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In the on-off model with nearest-neighbor hopping dynamics, the stationary 
distribution does not factorize and the stationary solution of the Master equation is 
not known. We therefore study the model using numerical Monte-Carlo simulations. 
From these simulations we find that condensation does indeed seem to be controlled 
by an effective hopping parameter b e g, albeit with a value which differs from the MF 
prediction. We also find that asymmetric jump rates may cause the condensate to drift 
with a finite velocity. 

In this section we concentrate solely on the on-off model with jump rates of the 
form ffl~4"l) and (125|) . unless explicitly stated otherwise. 

4-1. On-off model with symmetric nearest-neighbor hopping 

We begin the discussion of a ring with nearest-neighbor hopping dynamics by considering 
an on-off model with symmetric hopping, i.e., with p — 1/2. Note that, unlike the 
Markovian ZRP with symmetric hopping, which satisfies detailed balance and hence is 
an equilibrium model, the non-Markovian ZRP does not satisfy detailed balance even 
when it is symmetric. To understand why, note that there are allowed dynamical moves 
whose reverse cannot occur (such as an advance of a clock, or a jump of a particle 
simultaneously with resetting the clock of the target site to zero). As these moves have 
a non-zero probability to occur in the steady state, stationary probability currents must 
exist. 

Monte-Carlo simulations of the on-off model with nearest neighbor hopping were 
carried out on a ring of size L = 500 boxes with different particle densities and b = 4.5. 
The system was initialized to a state in which all particles were located at the first 
site and all sites were "on", and the dynamics was run for a time of i e quii = 2.5 • 10 7 
time units to allow the system to reach a steady state. After this equilibration time, 
the state of the system was recorded every Sampling = 5000 time units. The measured 
single site occupation probability P{n) and typical snapshots of the lattice for different 
values of p, presented in Fig. [2j show a qualitative resemblance to those of a Markovian 
ZRP. At small densities when the system is in the fluid phase, the single site occupation 
probability has an exponential tail, while at high densities this probability develops 
a peak which corresponds to the condensate. The transition takes place at a critical 
density (which is found to be p c ~ 2.75) at which the occupation probability decays as a 
power law of the form P{n) ~ n~ bcS . In finite systems this power law has an exponential 
cutoff due to finite size effects. 

Measuring the effective hopping parameter b e s numerically is a difficult task because 
it depends on the tail of the probability distribution which is strongly distorted by finite 
size effect. However, simulation results indicate that 6 e ff for symmetric hopping is larger 
than the MF value ( 1241) and smaller than the "bare" value b (see Figure [2]) . 

The simulation results indicate that the conclusions which were found for mean-field 
dynamics are qualitatively correct for symmetric nearest-neighbor dynamics. 



Motion of condensates in non-Markovian zero-range dynamics 



12 



10" 



°- 10" 



10" 



— Subcritical 

— Critical 

— Supercritical 



mean-field 



50 







Subcritical 




Uu li Li. ^lAJhjJduLkhJj 









100 



200 



300 



400 



500 



50 





Critical 



JLuttoauultJl ilUl ill J^jJLd^JiklLiM^JajLl. Jj 



bare 



10" 



10' 



10' 



10° 



1182 

60? 
30 




100 

Supercritical 



200 300 



400 



500 



JL i-Ui.iitu 



Hi,.. J, i 



lUilU. in 







100 



n 

(a) 



200 300 
Site /' 

(b) 



400 



500 



Figure 2: (a) The occupation probability P(n) of a single site in a ring with symmetric 
nearest-neighbor on-off dynamics for several densities, (b) Typical snapshots of the 
lattice for several densities. Results were obtained by Monte-Carlo simulations with 
b = 4.5, c = 1 and L = 500 sites, for a subcritical density (p = 2), a supercritical 
density (p = 5) and at the critical region (p = 2.75). The two thin straight line in 

of the bare value 



and P(n) 



~ n 



(a) correspond to the power laws P(n) ~ n~ 
b = 4.5 and the mean-field effective value 6 e fr ~ 2.78 of Eq. (I24j) . The power-law 
exponent for symmetric dynamics is seen to lie between these two values. Note that 
in the supercritical case, the y-axis of figure (b) is broken in order to show both the 
condensate and the disordered background. 



4-2. On-off model with asymmetric nearest-neighbor hopping 

Simulations of asymmetric nearest-neighbor dynamics (i.e., with p < 1/2) indicate 
that, as with the symmetric case, condensation is controlled by an effective hopping 
parameter. However, a new effect is found in simulations of asymmetric hopping: the 
condensate drifts with a finite velocity. Two different drift regimes are observed in 
simulations of finite systems: a "strong drift" regime in which the condensate is in a 
continual motion, and a "weak drift" regime in which the condensate stays for some 
(random) time in each site before jumping to the next. In what follows we begin 
by discussing the case of totally asymmetric hopping dynamics, (i.e., with asymmetry 
parameter p = 0), where we examine the strong drift and the weak drift regimes 
separately. We then discuss more general asymmetric dynamics, including partially 
asymmetric hopping and asymmetric dynamics on higher dimensional lattices. 

4-2.1. Totally asymmetric hopping: strong drift regime Monte-Carlo simulations of a 
ring with totally asymmetric nearest-neighbor on-off hopping dynamics were carried out 
for different values of c. The results of these simulations show that for small values of 
c the condensate drifts continuously in what we term a strong drift regime. In this 



Motion of condensates in non-Markovian zero-range dynamics 



13 




it 



6 



8 



— © t= 5-10 7 
— Q t= 10-10 7 
— t = 15-10 7 
— « t = 20-10 7 



□ 



4 



2 







35 



45 



245 



635 645 

site /' 




Figure 3: Snapshots of the on-off model with totally-asymmetric nearest-neighbor 
hopping on a ring, showing the occupation numbers n« at and in the vicinity of the 
condensate at four points in time. Here, L = 1000 sites, p — 10, b = 5.5, and c = 1. 
The condensate occupies two sites and drifts with a constant mean velocity. 

regime, the condensate typically occupies two adjacent boxes i and i + 1, in contrast 
to previously known condensation phenomena. In addition, The location of these two 
boxes advances with time. This is demonstrated in Fig. [31 where we present snapshots 
of the lattice taken at different times as obtained from a simulation with L = 1000 
boxes. 

An inspection of the microscopic dynamics shows that the drift of the condensate 
takes place via a "slinky" motion in which the second condensate site, i+l, accumulated 
particles at the expense of the first condensate site, i. This slinky motion results from 
the fact that site i + 1 is turned off more often than other sites. In other words, the 
effective hopping rates out of a site are no longer homogeneous in space, but rather 
they depend on the distance of the site from the condensate, and in particular, the 
mean current out of the condensate is larger than the mean current out of the next site: 
(ui(ni)) > (ui + x( n i+i)) ■ Thus particles accumulate on site i + l until site i is no longer 
macroscopically occupied, giving the clock at i + 1 the chance to reach the on state for 
durations of time sufficiently long to allow particles to escape. Then particles start to 
hop from site % + 1 to site i + 2 in the same fashion and the slinky motion continues. 
This mechanism for condensate motion was recently found in other models, and will be 
analyzed in more detail elsewhere [171 ] . 

This slinky motion mechanism suggests that the drift velocity thrift is inversely 
proportional to the number of particles in the condensate N cond , i.e., 



where iVbg = N — iV con d and pb g = ^bg/ (L — 2) are respectively the mean number and 
the density of particles in the background fluid, i.e., in all sites but the condensate sites. 
In the thermodynamic limit, the velocity of the condensate vanishes. It should be noted 
that this drift motion of the condensate is different from the relocation of the condensate 
which occurs in Markovian ZRPs. In the Markovian condensate on any finite 




(33) 
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Figure 4: (a) The "slinky" motion of the condensate is seen in the position of the most 
occupied site i max and its occupation number n max as a function of time, on a ring of 
size L = 1000 with c = 1. (b) This slinky motion suggests that the velocity of the 
condensate scales as L' 1 , as is demonstrated for different values of c. Graphs were 
obtained from simulation of totally-asymmetric nearest-neighbor on-off dynamics with 
p = 10 and b = 5.5. 



system can melt and reappear at some other randomly chosen site of the lattice. This 
relocation of the condensate happens on a characteristic time which scales with the 
system size to a power larger than 2 18N21|. A similar relocation of the condensate to 
a random distant site is seen to occur also in the asymmetric on-off ZRP, superimposed 
on the "slinky" drift motion. 

The snapshots presented in Fig. [3] clearly demonstrate that the condensate occupies 
two adjacent sites with varying relative occupation, consistent with the slinky motion 
described above. In addition, the drift of the condensate is evident in the figure. In 
order to demonstrate the slinky motion in more detail, we present in Fig. Ha] a plot 
showing the position of the most occupied site i max and its occupation number, n max , as 
a function of time. The occupation number n max oscillates in time with approximately 
constant frequency. Typically it decreases linearly until it reaches its minimal value, 
when i max increases by 1 and n max starts increasing. Fig. [4b] displays the scaling of the 
condensate velocity with the system size L, which agrees with the estimate of Eq. ( }33lh 

In Fig. [5] we present the single-site occupation probability distribution P{n) for 
various densities and for various system sizes. At high densities the distribution exhibits 
a plateau which reflects the particle distribution among the two sites which constitute 
the condensate. This is in contrast with a Markovian ZRP and the symmetric on-off 
model where the condensate is supported by a single site, which results in a sharp peak 
in P{n) (compare with the inset and with Fig. I2aj) . The value of P{n) at the plateau in 
the non-Markovian case may be estimated for p above the critical density and large L 
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Figure 5: The occupation probability P(n) of a single-site in a ring with totally- 
asymmetric nearest-neighbor hopping and different densities, as obtained from Monte- 
Carlo simulations, (a) P{n) at a subcritical density (p = 3), supercritical density 
(p = 10) and at the critical region (p = 4.1). The horizontal line indicates P p i a teau of 
( |34l) . where pb g was obtained from the simulation. Here, L = 1000, b = 5.5 and c = 1. 
The inset shows for comparison a similar plot of P(n) for a Markovian ZRP of L = 1000 
sites with b = 3, in the subcritical (p = 0.5) critical (p = 1) and supercritical (p = 4) 
phases, (b) P(n) for several system sizes (L increases in the direction of the arrow). 
The dashed horizontal lines indicates P p i a teau of ( l34"|) . Here p = 10, b = 5.5 and c = 1. 




using the slinky motion of the condensate. The probability that a given site carries the 
condensate is 2/L, and in such a site there is an approximately uniform probability to 
find any occupation < n < iV cond = L(p — p hg ). Thus, 

2 1 

-Pplateau ~ f T7 \ ■ (^4) 

L L{p - p bg ) 

This estimate is in good agreement with the plateau value in Fig. [5j For small densities, 
P{n) decays exponentially, indicating the absence of a condensate. For the system size 
studied in this figure, the distribution at small values of n does not allow to extract a 
power law decay as expected for the condensation transition. At density p = 4.1 there 
is a range of n for which P{n) seems to follow a power law with b e s ~ 4. This value 
differs significantly from the bare parameter b = 5.5, the expected value for Markovian 
ZRP. 

4-2.2. Totally asymmetric hopping: weak drift regime For larger values of the clock 
rate c, the motion of the condensate looks qualitatively different from that in the strong 
drift regime: the continuous slinky motion is replaced by an erratic slinky motion, in 
which the condensate spends a long period of time in each site before jumping to the 
next. We refer to this regime as the weak drift regime. This difference is observed 
on finite systems. Whether this type of motion persists for large L remains an open 
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(a) (b) 

Figure 6: (a) The position of the most occupied site z max and its occupation number 
Wmax as a function of time, in the weak drift regime, as measured on a ring of size 
L = 2000 with totally asymmetric hopping, b = 5.5, p = 10 and c = 2. The times Tdriftj 
T S piu and Tbamer are illustrated in the bottom panel, (b) The occupation probability 
P(n) of a single-site measured in the same system. A plateau is still observed, reflecting 
the slinky motion during transitions for one site to the next. However, unlike the strong 
drift regime, there is a peak at high densities reflecting the periods of time when the 
condensate is motionless. 

question at this point, as there are some indications that the thermodynamic behavior 
might be similar to the strong drift motion in this limit. In what follows we present the 
numerical evidence for the weak drift regime, and provide details on the question of the 
thermodynamic limit. 

All three main features which characterize the strong drift regime — a condensate 
that occupies two sites, its continual drift, and a plateau in the single site occupation 
probability — are modified in the weak drift regime. In this regime, the condensate 
occupies a single site for a long duration of time, and it occupies two sites only during 
the (relatively short) time of transition from one site to the next. This is clearly seen 
in Fig. [6aJ where the occupation and location of the most occupied site are shown as a 
function time for a system of size L = 2000 with c = 2 and b = 5.5 (compare with Fig. 
Hal and note the difference in the scale of the time axes). As a result, a sharp peak is 
seen in the single site occupation probability at large values of n (Fig. [6b|) . A plateau is 
still found at intermediate values of n, but it no longer follows the scaling relation ( 1341) . 

The characteristics of the weak drift described above suggest that the condensate is 
stabilized in one site by a "barrier" . The slinky motion is initiated only once fluctuations 
overcome this barrier and the number of particles in the condensate decreases beyond 
some threshold, or, alternatively, when the number of particles in the next site increases 
beyond a threshold. 

A possible microscopic mechanism which would give rise to such a threshold is as 
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Figure 7: (a) r bari . icr and r spill (see Fig. [Ha]) as a function of the system size, (b) The 
ratio Tspiii/Vbamer as a function of the system size. For the system sizes which we could 
simulate, the condensate motion has not yet converged to its thermodynamic limit 
behavior. If the trend shown here continues at larger L, eventual condensate motion 
will be similar to that seen in the strong drift regime. The parameters used in the 
simulation are the same as those of Fig. O 



follows. Suppose the condensate is located at site 1. As discussed above, the drift 
motion of the condensate indicates that the mean current out of site 1 is greater than 
that leaving site 2, i.e., (wi(ni)) > (u2(n 2 )). However, if n 2 is small enough, it might 
be that at some moment (wi(ni)) < U2(n 2 ). In this case, because Ui{n) is a decreasing 
function of n, particles begin to accumulate in site 2 only after its occupation exceeds 
a value n* which is defined by v,2{n*) = (ui(ni)). Thus, the condensate begins to spill 
from site 1 to 2 only after a random fluctuation brings the occupation of site 2 to n*. If 
n* is large enough, the time until such a fluctuation occurs can be long. However, this 
time is expected to remain finite in the thermodynamic limit L — > oo. If this picture is 
correct, the erratic motion of the condensate which characterizes the weak drift regime 
is expected to be negligible in the thermodynamic limit, since the time of the spilling of 
the condensate scales as the system size L. A more detailed study of such a mechanism 
for a weak condensate drift will be presented elsewhere (l7| . 

It is not yet known whether this picture provides an accurate description of the 
microscopic mechanism which leads to the weak drift motion. However, numerical 
evidence indicates that the weak drift regime may indeed exist only as a finite size effect. 
To address this question, we compare the typical time that the condensate resides on a 
single site, which we term Tbarrien with the time it takes the condensate to "spill" from 
one site to the next, which we denote r sp iu- Together, these two time add up to give 
the typical time for the drift motion: r drift = t>^ r - ft = r spill + r barrier , see Fig. [6a] In Fig. 
ITal we present the dependence of rbarrier and r sp m on the system size. For the system 
sizes which we were able to study numerically, r sp in was seen to grow linearly with L as 
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Figure 8: The mean time between (a) drift movements of the condensate (rdnft) and (b) 
melting movements (r me it) as a function of c for several values of L. The entire duration 
of the simulation was ~ 10 9 time units, and thus the plots are cut off at r ~ 10 8 , beyond 
which the melting time can no longer be estimated reliably with the data available. The 
parameters of the simulation are the same as those of Fig. [6l Lines are guides to the 
eye. 



expected (it should takes twice as long to move twice as many particles from one site 
to the next). However, Tbamer is seen to grow slower than linearly. This trend, which is 
emphasized when looking at the ratio r spi u/r barriei . (see Fig. l7b~|) indicates that although 
Tspin <C Tbarricr for the system sizes which were studied, the situation might be reversed 
at large enough systems, in which case the motion of the condensate will be similar 
to that in the strong drift regime. Whether this trend continues at larger values of L 
remains an open question. 

Numerical limitations also hindered the study of the behavior of the system at the 
transition between the strong and weak drift regimes, as well as at higher values of c. 
At values of 1 < c < 2, there is a sharp decrease of the mean time r me i t between events 
at which the condensate melts and reappears in a distant site (see Fig. l8al) . For the 
values of L which we were able to study these events were still quite frequent, indicating 
that the system was still far from thermodynamic behavior. As Thrift is seen to grow 
roughly exponentially with c (see Fig. l8bl . when c is larger than about 2, rdnft becomes 
comparable with the total length of the simulation. 

4-2.3. Other types of asymmetric dynamics The main features of the on-off model, and 
specifically the drift of the condensate which was discussed above for the case of totally- 
asymmetric hopping, are quite robust to small changes in the dynamics of the model. 
We shall now mention a few such modified models which exhibit a similar behavior in 
the condensed phase. 

We begin with the on-off model with partially asymmetric dynamics, where each 
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Figure 9: Drift of the condensate in a 2- dimensional square lattice with bias to the 
up and right directions. The top two panels show the x- and y-components x max and 
Z/max of the location of the most occupied site as a function of time and the bottom 
panel its occupation number n max . The lattice is of size 20 x 20, and the parameters of 
the simulation are p = 30, b = 5.5 and c = 0.4. The asymmetry of the hopping is as 
described in the text. 



time a particle hops it can jump to the right with probability 1 — p or to the left with 
probability p. Totally asymmetric dynamics corresponds to p = 0. If an asymmetric 
system is in the strong drift regime and p is increased slightly, no significant changes 
in its behavior are seen, and in particular it remains in the strong drift regime. When 
p is further increased, a transition to the weak drift regime occurs in the numerical 
simulations. This transition is similar to the one discussed above in the totally- 
asymmetric case when c is increased beyond 1, and it too is accompanied by a sharp 
dip in r me it. For a system of size L = 1000 with c = 1 and b = 5.5 the transition was 
found to occur at around p = 0.1. Beyond this transition, the drift velocity rapidly 
decreases as the dynamics approaches the symmetric dynamics at p = 1/2 at which 
point no drift of the condensate is seen. It should be noted that in the symmetric case, 
when the condensate relocates to a different site there seems to be no preference to its 
neighboring sites. Rather, the condensate melts and reappears at a distant site, as in 
the Markovian case. 

A drift of the condensate is also observed when the site is not turned completely 
off at t = 0. Simulations with hopping rates of the form u{n,r) = u{ri)v{r) with v(t) 
as in Eq. (l32l) . vq = 0.5 and b = 5.5, exhibit strong drift behavior when c = 0.4 and 
weak drift behavior when c = 1. 

Condensate drift, both weak and strong, also occurs in 2-dimensional nearest- 
neighbor asymmetric on-off models. A particularly interesting case is when the hopping 
bias is not parallel to any of the lattice directions. Fig. [9] displays the motion of the 
condensate in a 20 x 20 square lattice with periodic boundary conditions where each 
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time a particle hops it either moves one site up or one site to the right with equal 
probabilities. The figure shows the x- and y-coordinates of the most occupied site and 
its occupation. It is intriguing to notice that the condensate moves alternatively up 
and to the right in quite an orderly fashion. Snapshots of the lattice (not presented 
here) reveal that the condensate typically consists of an L-shaped group of three highly 
occupied sites. At higher values of c the orderly motion is destroyed and the condensate 
drifts in the weak regime. 

5. Exactly solvable non-Markovian ZRP 

In this section we present a non-Markovian ZRP whose steady-state probability 
distribution factorizes into single site terms, similar to the usual Markovian ZRP. Thus, 
the steady-state distribution, the effective hopping rates and 6 e ff can be calculated 
exactly. A version of the model with totally asymmetric hopping is analyzed first in 
Sec. 15. 1^ and then the model and results are generalized to symmetric and partially 
asymmetric hopping in Sec. 15.21 

5.1. Totally asymmetric dynamics in 1-d 

5.1.1. Description of the model The exactly solvable non-Markovian ZRP is a variant 
of the on-off model in which the advance of a clock of a site depends on the clock 
states of neighboring sites. Before presenting general results for partially-asymmetric 
hopping and for lattices in any dimensions, we begin for simplicity by considering a 
one-dimensional lattice with totally-asymmetric hopping. 

The model is similar to the one described above in Sec. [2j at each site i of a 
one- dimensional lattice of L sites there are Hi particles, and a clock variable T\ = 0, 1, 
signifying "on" and "off". Note that our notation here differs from that of Sec. 13. 2\ 
where r was allowed to take any integer value. This is not a significant difference, since, 
as discussed there, identifying all r > 1 clock states with r = 1 does not affect the 
dynamics. A particle can hop from site % to % + 1 with rates ( fl4l) . and once a particle 
jumps the clock at the target site is reset to zero. The only difference in the dynamics of 
the exactly solvable model is in the way the clocks are updated: the clock at site i can 
change from to 1 only if Tj+i = 1. The allowed dynamical moves can be summarized 
as 

...,(ni,l),(n i+1 ,r i+ i),... U -^4 . . . , (rii-1, 1), (n i+ i + l,0), . . . 
...,(ni,0),(n m ,l),... ...,(rii,l),(n i+1 ,l),... (35) 

(here (rij,Tj) signifies the occupation and clock state of site i). 

5.1.2. The steady-state distribution The goal of this section is to construct the steady 
state distribution of this model and show that it has a factorized form. Before doing so, 
we note that the factorized form is somewhat different from that presented in Eq. (j3J). 
The reason for the difference is that states in which all sites are off cannot be reached 
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by the dynamics of the model (all other states are possible). This introduces some 
correlations between the sites beyond those generated by the conservation of particles. 
The product measure which we discuss below is therefore of the form 

L L 
i=l i=l 

where S(t) = 1 if r = 0, i.e., if all r« = 0, and is zero otherwise. Here f(n,r) are the 
single-site weights. The normalization is accordingly given by 

L L 
n,T i=l i=l 

In the thermodynamic limit, the weight of configurations with r = becomes negligible, 
and therefore adding the square-brackets term in fl36l) and fl37|) does not affect this limit. 
We describe some properties of this factorized form in Appendix B Note that the same 
product form also describes the finite-size product measure of the mean-field on-off 
model which was considered in Sec. 13.21 

The dynamics fl35|) defines an ergodic process (on the set of all configurations with 
a given number of particles and at least one "on" site), and therefore it has a unique 
steady state distribution. We now show that this distribution has a factorized form 
(|36|) . This is done by assuming such a factorized form, and showing that it is indeed 
the unique stationary solution of the master equation. We begin by writing down the 
master equation. To this end we define a function W ruTi+1 (ni,n i+ i) by 

V(n, r)Wi,i(ni, n i+1 ) = V(n, {. . . , n-i, 0, 1, . . .})c - V(n, r)w(nj), 
V{n,T)Wo,o{ni,n i+1 ) = 0, 

V(n, T)W , x {n h n i+1 ) = - V(n, r)c, (38) 
V{n,T)W 1>0 {ni,n i+1 ) = - V{n, T)u{rii) + ufc + 1) x 



P[{. 



r'=0,l 

, Hi + 1, n i+ i -!,...},{..., Tj_i, 1, r', r i+2 , . . .} 



Note that W rijTi+1 (rij, n i+1 ) in fact depends on the full configuration (n, r). We suppress 
this dependence in the notation because if V(n, r) has a factorized form, W indeed 
depends only on the occupation and clock states of two adjacent sites (see Eq. (TJ4]) 
below) . 

Using the function W, the master equation can be written as 

L 

P(n,T)=P(n,T)]TwV i+1 ( (39) 
i=i 

We elucidate Eqs. ( 1381) and ( |39i) through an example. Consider a configuration of a 
lattice of 4 sites with clocks r = {1,1,0,0} and some occupations n. The different 
transitions into this configuration and out of this configuration can be enumerated one 
bond at a time: 
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• Sites 1 and 2: Since both clocks are on, the only possible transition involving both 
sites which would lead to this configuration is an advance of the clock at site 1, 
which occurs with rate c. The only possible transition out of this configuration 
which involves the two sites is a particle hopping from 1 to 2. Therefore, this bond 
contributes two terms to the master equation, 

V(n, {0, r 2 , . . .})c - V(n, r)u(m) = V{n, T)W 1A {n u n 2 ). (40) 

• Sites 2 and 3: The first clock of the two is on while the second is off. The 
only possible transition involving these two sites leading to this configuration is 
a particle hopping from 2 to 3, and this is also the only possible transition out of 
this configuration. Therefore, this bond contributes to the master equation 

£ V({...,n 2 + l,n 3 -l,. ..},{...,!, r 3 , . . .})u(n 2 + 1) - V{n, r)u{n 2 ) = 

T 3 =0,1 

= V(n,T)W lt0 (n 2 ,n 3 ). (41) 

• Sites 3 and 4: Both clocks are off, and therefore no transitions which involve only 
this bond are possible. One can define the contribution to the master equation as 

= V(n,r)W 0fi (n 3 ,n 4 ). (42) 

• Sites 4 and 1: The first clock is off and the second is on. There are no transitions 
involving only this bond which can lead to this configuration. However, there is a 
possible transition out of this configuration, by an advance of the clock of site 4. 
The contribution from this bond is therefore 

-V(n,T)c = V(n,T)W ,i(n±,n 1 ). (43) 

Summing up Eqs. fj40l - (f43|) leads to the master equation (l39l) . A similar analysis shows 
that the master equation has exactly the same form for any configuration and for any 
lattice size. 

Now assume that P(n, r) has the factorized form fl36l) . In this case, the definition 
( |38l) has the simpler form 

W ltl (ni, n i+1 ) = c - u{ni), 

W Q> o(ni,n i+1 ) = 0, 

W 0il (ni,n i+1 ) = - c, 

ur i \ fon(ni + l)f(n i+ i-l) 

Wi fl {n h n i+1 ) = u{rii + 1) - u(ni), (44) 

jon{ni)jo&{n i+ i) 

where f oS (n) = /(n,0), f on (n) = f(n, 1), and f(n) = E T '/( n ' r ')- In the 
steady state, the left-hand side of Eq. (|39l) vanishes and the equation becomes 
J2i=i W Ti ,T i+1 {nii n i+i) = 0. This equation is solved by explicitly constructing its unique 
solution. This is done in two steps. First, we show that if one finds f(n,r) which 
satisfies 

= ^1,1(714,^+1) (45) 
= Wi t0 {n h m+i) + ^0,1(71^,71^+1) (46) 
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for any n i} n i+ i, rij, rij + i, this / is a solution to the equation. Then, we construct such 
an /. 

The first step is achieved by noting that the number of Wi t o(rii,ni+i) terms in 
the sum (139]) exactly equals the number of Wq^Uj, rij+i) terms in the sum, since any 
configuration of r^'s must have the same number of 01 and 10 nearest-neighbor pairs. 
Therefore, all terms in the sum (139~1) vanish either individually or in pairs, and the sum 
equals zero. We now construct a solution / which satisfies (j4"5]) - (j4"S]) . Condition fj45|) is 
equivalent to 

cfos(n) = u(n)f on (n). (47) 

For condition (1461) . note that Wo,i( n i) n j'+i) * s m f &c t independent of Uj,nj + i (see Eq. 
(144j) ). Therefore, this condition together with (1471) yield 

HjH + 1) u(nj + 1) = f(n i+ i) u(n i+1 ) 
f(ni) c + u(rii + l) f(n i+ i - 1) c + u(n i+ i) ' 
As the occupations rii and n^+i may vary independently, this equation holds only if both 
sided are equal to a constant, which might be set to 1/c without loss of generality (as 
it only affects the normalization constant Z^n). We therefore find that 

f(n)u(n) = f{n-l) (49) 

where 

111 

TT = - + TV ( 5 °) 
u[n) c u[n) 

compare with equations (jUJ) and (|2D|) . Choosing the constant to be 1/c guarantees that, 
as we show below, u(n) as defined in Eq. (150]) are the effective hopping rates. 

The conclusion from Eqs. (l4Tj) - (!50l is that the factorized probability distribution 
of the form (I3(jp with 

f{n, 1) = f on (n) = J . , f(n), (51) 
c + u[n) 

/(„,0) = / offi (n) = -^ T /(n), (52) 
c + u[n) 

n ., n 1 1 

/(n) =n^=n[- + ^]- («) 

is the stationary solution of model. It is easy to verify using (IBT?]) and (15T]) that 
u(n) = ^ on ^"( n ) ; and therefore u(n) are the effective hopping rates as defined in Eq. 



{TQj). 

Using the results ( |5T]) -( l53|) . one can calculate numerically the stationary probability 
for any configuration in a finite system of size L with N particles (recursion relations 



that facilitate this calculation are presented in |Appendix B[ ). Condensation in the model 
is determined by the probability measure in the thermodynamic limit. The factorized 
product measure ( |36i) and Eq. ( 1491) imply that this model has the same thermodynamic 



behavior as a Markovian ZRP with effective hopping rates (|50|) (see Appendix B ). One 
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can thus study condensation in the model using known properties of the ZRP, as was 
done in Sec. [3j In particular, for rates of the form u(n) = 1 + b/n + 0(n~ 2 ), one finds 



c l + £ + Q(n- 2 ) 
uin) = £ — — = J c 

1 + c l+(T^ + 0(^ 2 ) 



VOin ) 

n 



(54) 



with 



J c = < 1 and b e s = J c b < b. (55) 

1 + c 

Here, J c is the current at the critical density, and 6 e ff is the parameter controlling 
condensation. In other words, condensation may occur only when fe eff > 2, or b > 2/J c 
(compare with the MF values in Eqs. (I21~]) - (l24")) . and note also that at criticality, the 
probability to find a site in the off state is P Q s = J c /c, rather than ( flTl) of the MF 
model). 

5.2. Partially asymmetric dynamics and higher- dimensional lattices 

5.2.1. Description of the model The exactly solvable model described above can 
still be fully analyzed when the dynamics is generalized to partially asymmetric or 
symmetric dynamics and to certain higher dimensional lattices. Moreover, the stationary 
distribution turns out the be independent on the asymmetry or the dimension. We 
described the generalized dynamics and its solution in this section. 

First, consider dynamics on a 1-d lattice that allows for partially asymmetric 
hopping. This is implemented as discussed above in Eq. fl3]): when a particle jumps 
from an "on" site i (an event which occurs with a rate u(n)), it randomly chooses 
its target site: with probability 1 — p it moves to site i + 1 and otherwise (i.e. with 
probability p) it moves to % — 1. Here < p < 1 is the asymmetry parameter: p = 1/2 
corresponds to symmetric dynamics, while p = corresponds to a totally asymmetric 
bias to the right. 

For the stationary distribution to factorize, one must also modify the update rule 
for the clock variable, in the following manner. At each "off" site, an attempt to update 
the clock is made with rate c. Once an attempt is made at, say, site i, a neighboring site 
is chosen at random with the same asymmetry parameter p: site i + 1 is chosen with 
probability 1 —p and site i — 1 with probability p. Finally, if the chosen site is in an "on" 
state, the clock of site i is turned on. The generalized dynamics can be summarized as 

. . . , (m, 1), (rij+i, n + i), ... (1 ^H ni) . . . , (rii-1, 1), (n i+ i+l, 0), . . . 

. . . , (nj_i,Tj_i), (n^ 1), • • • ■ ■ ■ , (rji_i + l,0), 1), . . . 

...,(ni,0),(n i+1 ,l),... ...,(m,l),{n i+1 ,l),... 

...,(rai_i,l),(ni,0),... ^> . . . , (ni_i, 1), (n^ 1), • • ■ , (56) 

where the first line describes a particle jump from i to the right, the second describes 
a jump to the left, and the third and fourth lines describe the two update processes of 
the clock at site i. 
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In a similar fashion, the model can be generalized to symmetric or biased dynamics 
on higher dimensional lattices. Here we consider for concreteness cubic lattices in d- 
dimensions, although the argument which we present below for the factorization of 
the stationary distribution is valid for other lattices, e.g. a triangular lattice in 
As in the partially asymmetric particle leaves any site i, if it is on, with rate 

u(n). It then selects its target from among the 2d nearest neighbors of i according 
to an asymmetry probability vector p a , where a = 1, . . . , 2d denotes the direction (for 
example, in two dimensions a = 1,2,3,4 could correspond to north, east, south and 
west) and Y^=iPa = 1. A choice of p a = l/2d for all a corresponds to symmetric 
dynamics, and any other choice would result in biased hopping. The clock update rule 
in the d- dimensional case is similarly generalized: if site i is off, an attempt to update 
its clock is made with rate c. At each attempt, the neighbor of i in the direction a is 
chosen with probability p a , and if the chosen neighbor is on the clock of i is updated. 

Since 1-dimensional partially asymmetric hopping is a particular case of d- 
dimensional dynamics, both cases are treated below together. 



5.2.2. The steady-state distribution The factorization in the generalized case is 
demonstrated as done above, by explicitly constructing the stationary measure. To 
this end it is once again assumed that the stationary measure has the factorized form 
( |36]) - (I3T1) . The master equation for this factorized distribution reads, at the steady 
state, 

2d L 

= V(n, r) = V(n, r) W n , n+a {n h n l+a )\ , (57) 

a=l i=l 

where site i + a denotes the neighbor of site i in the direction a, and W is defined in 



The key observation which facilitates finding a solution to this master equation is 
that for each i and a such that T{ = 1 and r i+a = 0, there exists exactly one site j whose 
clock is Tj = while Tj +a = 1. This can be seen for example by examining the clocks 
of all sites on the ray which starts at site i and is in direction a (i.e., by examining 
sites i + 2a, i + 3a, . . .), which leads to a situation similar to the one- dimensional case. 
Therefore, a solution to Eq. (1571) can be found if Wi,i(n, n') = and Wi$ (n, n')+Wo,i = 
for all n and n' (note again that W^i = Wo t i(n,n') is independent of n,n'). These are 
precisely the conditions which appeared in the totally asymmetric case, and therefore 
they are fulfilled by the same solution — Eqs. (l5Tl) - (!53l) . 

We have thus shown that the stationary distribution of the generalized model 
factorizes, and moreover it is independent of asymmetry and lattice dimension. In 
particular, condensation is independent of the asymmetry parameter, and the results of 
Sec. 15.1.21 apply. 



| Note, however, that the argument does not hold for all higher dimensional lattices. For example, the 
argument fails for a 2d honeycomb lattice. 
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6. Conclusions 

The analysis presented above reveals that non-Markovian dynamics may have two major 
effects on the condensation transition of the ZRP. First, the parameter b which controls 
condensation is "renormalized" by the existence of memory in the dynamics, and thus 
a memory may suppress or induce condensation. For models with mean-field dynamics 
and for an exactly solvable variant of the model, the effective rates could be computed 
exactly, and thus the modified criterion for condensation was found. Numerically, the 
condensation in models with nearest-neighbor hopping were also found to be controlled 
by an effective b, although one which differs from the mean-field value. Calculating 
the effective hopping rates in nearest-neighbor models remains an open problem which 
may be of practical importance when one wishes to use a non-Markovian ZRP to study 
condensation in other systems. 

A second effect of the memory is perhaps more dramatic: the condensate is found 
to move from one site to the next when the dynamics is of asymmetric nearest-neighbor 
hopping. Numerical studies of finite systems identify two modes of condensate drift: a 
strong-drift regime with continuous "slinky" motion and a weak-drift regime in which 
the motion is more erratic. Both modes of motion are rather robust to changes in the 
dynamics. The behavior of the model in the thermodynamic limit is not yet known, and 
it would be interesting to ascertain whether there is a sharp transition between them, 
or, if such a transition does not exist, to understand the crossover from one regime to 
the other. 

The mechanism which leads to the condensate drift is understood on a heuristic level 
and is expected to be a generic feature of many systems which undergo a condensation 
transition and which are asymmetric and have some spatial correlations jljj]. However, 
a more quantitative understanding of this drift, for example the calculation of the drift 
velocity, remains an important open problem. It is also interesting to explore similar 
effects in other mass-transport systems, such as driven diffusive systems and shaken 
granular gases. In this respect, it should be noted that a mass-transport model with 



a moving condensate was recently identified in 22]. There, a variant of the ZRP is 
studied which has a factorized steady-state and in which unbound hopping rates lead to 
a condensate which reaches an infinite velocity. A product measure steady state with a 
moving condensate is not possible in systems with finite hopping rates like ours. 

We have also studied an exactly soluble variant of the non-Markovian model with 
nearest-neighbor hopping whose steady state factorizes. In this variant, as in the mean- 
field model, condensation is controlled by an effective b and no condensate motion 
appears. It should be notes that although the model has a product measure, particle 
currents are temporally correlated. 
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Appendix A. Regular and irregular clocks 

As stressed above, the internal clock variables Tj do not measure an exact time, but 
rather proceed in an irregular stochastic fashion. In this Appendix, it is shown that 
regular clock, that proceed in a deterministic continuous fashion, may be obtained from 
the dynamical rules (PQ) by taking an appropriate limit. 

We denote the clock variables in this Appendix as m; instead of r^, to emphasized 
that they may attain only integer values. In order to obtain regular clocks, define new 
clock variables T{ = m^dr, where dr is an infinitesimal time unit which will eventually 
be taken to zero. The new clock variables are no longer integer: they can attain any 
value r = 0, dr, 2dr, 3dr, . . ., and in the limit of infinitesimal dr they become continuous 
variables. In addition, the rate with which m; advances to + 1 is taken as c = 1/dr. 
Finally, the hopping rates out of each site i are taken to depend on Tj rather than mi, 
and thus they can be written as u(ni,Ti). The limit of regular clocks is then obtained 
by taking the limit dr — > 0, while keeping Tj fixed. 

For example, consider an on-off model with regular clocks, whose hopping rates 
are u(n,r) = u(n)Q(r/To — 1), where To is a constant and Q(x) is the Heaviside theta 
function. Such a model may be achieved by considering irregular-clock models (CQ) with 
rates u(n,m) = u{n)Q{m dr/T — 1) and c = 1/dr, and taking the limit dr — > while 
keeping the constant To fixed. In this regular-clock on-off model, whenever a particle 
hops into a site this site is turned off for a duration of exactly r time units. The 
solution of such a model with mean-field dynamics may be found from an analysis 



similar to that presented in Sec. |3] 23]. Similarly, more general regular-clock models 
with rates u(n, r) = u(n)v(r) can be obtained by taking the limit of irregular-clock 
models with rates u(n,m) = u(n)v(mdT /tq) (the function v remains unchanged when 
taking the limit). 



Appendix B. Properties of factorized distributions of the form (136ft 

In this Appendix we present some of the properties of the stationary distribution of 
the exactly solvable on-off model, which has the factorized form ( )36|) with partition 
function ( 137]) . The goal of this Appendix is to present recursion relations which allow 
the calculation of this product measure for any finite system size, and to demonstrate 
that such product measures lead to the same thermodynamic behavior as fll]) and ([5]). 
We begin by defining two auxiliary partition sums, 

L L 

n,x i=l i=l 
L L 

z °l s n= En^K5> _iv )' (R2) 

n i=l i=l 
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and two auxiliary distributions, 

L L 

P(n, r) = _ TT f(rn, n) <*(V n, - iv), (B.3) 

L L 

V° S (n) =^ r \{h{n l )b-(^n l -N). (B.4) 

L,N i=i i=i 

Here and in the rest of this section we denote fo(n) = f(n,0) and /i(n) = f(n,l). 
This is done to avoid confusion with the superscript "off", which will be used below 
to denote quantities calculated using the distribution (1B.4jl . As before, we denote 
f(n) = fo(n) + fi(n), and we adopt the convention of Sec. [5] whereby the clocks may 
have only two values, r = 0, 1. 

Using these notations and the definition (137)) one immediately finds that 

Zl,n = Z L:N — Z°^ N . (B-5) 

We first analyze the auxiliary distributions before treating the original problem. By 
summing over the occupations and clock states of all sites but one, the probability to 
find a single site in any given state is found to be 

L 

P(n,r)= ^2 n 2 , . . . , n L ; r, r 2 , . . . , t l )S (^] m - (N - n) 



n 2 ,...,n L j_ 2 



/(n,r)%^, (B.6) 



P oS (n) = Yl -P° S (n,n 2 ,...,n L )6(22ni-(N-n))=f (n) . (B.7) 

n 2 ,...,n L i=2 L > N 

Summing both equations over n leads to the recursion relations 

N 

Zl,n= y^L-^N-nfjn), (B.8) 



n=0 
N 

70S _ 
J L,N — 

n=0 



E^V-n/o(n). (B.9) 



When f(n, r) are known, these recursion formulas can be used for a numerical calculation 
of the auxiliary partition sums, and thus, using (1B.5j) also of Zl,n- 

Knowing the partition function Z^n, other quantities of interest can be computed. 
For example, repeating the calculation of (1B.6)) for the original product measure yields 

TDl \ ft \ Z L-l,N-n + br,lZfL lN _ n 

P{n, t) = fin, t) ! , (B.10) 

where S Tt \ is the Kronecker delta. The current can be found in a similar fashion by 
calculating 

N N - N—l - 

(«(n)> = P ^ !)«(») = E = E m Zh -y N ~ l ~ n 
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= (B.H) 

where fl5Tl) - fl53|) were used to deduce that f(n,l)u(n) = f(n — 1), and we have used 
f lR5|) and flRlQ]) . 

In the thermodynamic limit, the partition function can be analyzed by transforming 
to the grand-canonical ensemble. Mathematically this is done by introducing the grand- 
canonical partition function which is the generating function 

oo 

Z L {z) = ^2z N Z LtN . (B.12) 

N=0 

Using the definition (1371) . one can split the sum into two contributions, Zl(z) = 
Zl — Zf 1 , where 

oo oo 

Z L = z N Zl,n = F(z) L , Zf = ^Kn = F oS {z)\ (B.13) 

iV=0 N=0 

and 

oo oo 

F(z) = f(n)z n , F° s (z) = fo(n)z n . (B.14) 

n=0 n=0 

Using (15ip — (155)1 . one has F oS (z) = zF(z)/c, from which the grand-canonical partition 
function is found to be 



Zr(z) = Zt(z) 



1 

c 



(B.15) 

If the radius of convergence of the sum (IB. 12)) . or equivalently of (IB.14I) . is smaller than 
c, then the correction due to the weak correlation between clocks is exponentially small 
when L is large. For hopping rates of the form u(n) ~ 1 + b/n, this radius of convergence 
is z c = c/(l + c) < c (see Eq. (153)) ). and therefore, in this case (z/c) L is indeed negligible. 
Note that z c is the current J c of the canonical system at the condensation transition, 
Eq. ([55]). 

The relation between the fugacity z and the canonical density p is given by 
* dlog Z(z) F'(z) 1 

p= i^z— =z m (R16) 

which is an implicit equation for z(p). This is the same expression as that of a Markovian 
ZRP with rates u(n), up to a correction which is exponentially small in L. 
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